Device and method for multispectral optoacoustic imaging

ABSTRACT

The invention relates to a device and an according method for multispectral optoacoustic imaging of an object, the device comprising: an irradiation unit configured to irradiate the object with electromagnetic radiation at two or more different irradiation wavelengths (λ), said electromagnetic radiation having a time-varying intensity; a detection unit configured to detect acoustic waves generated in the object upon irradiating the object with the electromagnetic radiation at the different irradiation wavelengths (λ); and a processing unit configured to reconstruct images of the object based on the detected acoustic waves generated in the object at each of the irradiation wavelengths (λ) and to determine a spatial distribution of at least one first concentration value, which relates to a concentration of at least one electromagnetic radiation absorber in the object, wherein the determination of the spatial distribution of the at least one first concentration value is based on the reconstructed images at the different irradiation wavelengths (λ), on at least one wavelength-dependent extinction coefficient of the at least one electromagnetic radiation absorber in the object, and on a linear combination of at least two model spectra representing wavelength-dependent basis functions of a radiation fluence or of a function of a radiation fluence, in particular a normalized radiation fluence, in the object.

The present invention relates to a device and a method for multispectraloptoacoustic imaging, in particular for determining a quantitative imagerelating to concentrations of one or more chromophores in an object,according to the independent claims.

Optoacoustic imaging is based on the photoacoustic effect, according towhich ultrasonic waves are generated due to absorption ofelectromagnetic radiation by an object, e.g. a biological tissue, and asubsequent thermoelastic expansion of the object. Excitation radiation,for example laser light or radiofrequency radiation, can be pulsedradiation with short pulse durations or continuous radiation withvarying amplitude or frequency.

An important advancement in optoacoustics is imaging of differentmoieties by resolving their distinct absorption spectra. This istypically achieved by illuminating the object imaged at multiplewavelengths and performing spectral unmixing operations in the collecteddata. Imaging of such moieties is important to different applications.

In biological and medical applications for example, tissue bloodoxygenation is a significant metric that can infer valuable knowledgeregarding tissue viability, metabolism, hypoxia and neuronal activation.Few imaging modalities can provide quantitative information of deeptissue blood oxygenation in high spatial and temporal resolution. Byutilizing several excitation wavelengths so-called multispectraloptoacoustic images (also referred to as spectroscopic optoacousticimages) of high spatial and temporal resolution can be obtained, thatare related to tissue optical absorption, which can be used to inferblood oxygenation. However, accurate quantification of a spatialdistribution of chromophores, like oxygenated and deoxygenatedhemoglobin, in deep tissue by means of multispectral optoacoustictomography (MSOT) remains a challenge, in particular at depths exceeding3 mm.

Optoacoustic images are proportional to the absorbed energy densityH(r,λ) (where r denotes the position and λ denotes the wavelength of theexcitation source), which depends on the spatial distribution of theoptical absorption coefficient μ_(α)(r,λ) of tissue and the opticalfluence φ(r,λ), i.e. H(r,λ)=μ_(α)(r,λ) φ(r,λ). Assuming μ_(α)(r,λ)known, the concentrations of tissue absorbers can be estimated in eachposition r using linear spectral unmixing methods. Thus valuablephysiological (e.g. blood oxygenation in tissue) and molecularinformation (e.g. spatial localization and concentration estimation ofmolecular probes) can be extracted. However, the μ_(α)(r,λ) distributionis not explicitly known from raw optoacoustic images, since theinformation obtained is coupled to the unknown optical fluence φ(r,λ).Decomposing the absorbed energy density H(r,λ) into the optical fluenceφ(r,λ) and the optical absorption component μ_(α)(r,λ) is a majorchallenge in the MSOT problem.

A possible formulation for solving the problem of quantitativeoptoacoustic imaging (also referred to a quantitative photoacoustictomography, or qPAT) is a model-based optical property estimationapproach. Under this formulation, the optical fluence φ(r,λ; μ_(α)(r,λ),μ_(s)(r,λ)) is related to the spatial distribution of the opticalabsorption μ_(α)(r,λ) and scattering μ_(s)(r,λ) coefficients of thetissue imaged, using a light propagation model. A number of non-linearmodel-based inversion methods have been published that seek tosimultaneously extract the optical fluence and tissue optical properties(absorption and scattering coefficients), given the absorbed energydensity H(r,λ) at a single or multiple wavelengths by minimizing thefunctional ∥H(r,λ)−φ(r,λ; μ_(a)(r,λ), μ_(s)(r,λ)) μ_(α)(r,λ)∥. Opticalproperty estimation methods may operate in the signal or the imagedomain, they may operate under one excitation wavelength (where opticalscattering needs to be known for achieving uniqueness) or under multipleexcitation wavelengths and they may also be combined with a model foracoustic propagation. Their common characteristic is that they seek tominimize the functional ∥H(r,λ)−φ(r,λ; μ_(α)(r,λ),μ_(s)(r,λ))μ_(α)(r,λ)∥ for obtaining quantitative values of tissueabsorption coefficients.

A major disadvantage of the optical property estimation approach is thatit typically depends on the accurate estimation of the absorbed energydensity H(r,λ) in the whole illuminated tissue volume. Typically,optoacoustic images P(r,λ) are proportional to the absorbed energydensity H(r,λ) and to a space dependent scaling factor C(r), i.e.P(r,λ)=C(r)H(r,λ), where C(r) depends on the spatial impulse response ofthe system and the (generally unknown) spatially varying Grüneisencoefficient. Thus, the estimation of H(r,λ) from optoacoustic images maybe prone to error due to scaling effects typically not included intomographic reconstruction (such as the spatial impulse response of thesystem or ultrasound attenuation) or uncertainties in estimating theGrüneisen coefficient of different tissue types. Moreover, if theabsorbed energy density is unknown or significantly wrong in one part ofthe illuminated volume (common in the case of limited view systems orunder ultrasound scattering incidents) the inversion process may becompromised. Finally, optical property estimation approaches definelarge scale inverse problems with thousands of unknown parameters, whichresult in high computational complexity. For these reasons theapplication of such methodologies in the case of experimental tissuemeasurements of complex structure, such as biological tissues, is so farlimited.

It is an object of the invention to provide a device and method formultispectral optoacoustic imaging allowing for an improvedquantification of a spatial distribution of at least one value relatingto the concentration of one or more chromophores in an object.

This object is achieved by the device and the method according to theindependent claims.

The device for multispectral optoacoustic imaging of an object accordingto the invention comprises: an irradiation unit configured to irradiatethe object with electromagnetic radiation at two or more differentirradiation wavelengths, said electromagnetic radiation having atime-varying intensity; a detection unit configured to detect acousticwaves generated in the object upon irradiating the object with theelectromagnetic radiation at the different irradiation wavelengths; aprocessing unit configured to reconstruct images of the object based onthe detected acoustic waves generated in the object at each of theirradiation wavelengths and to determine a spatial distribution of atleast one first concentration value, which relates to a concentration ofat least one electromagnetic radiation absorber in the object, whereinthe determination of the spatial distribution of the at least one firstconcentration value is based on the reconstructed images at thedifferent irradiation wavelengths, on at least one wavelength-dependentextinction coefficient of the at least one electromagnetic radiationabsorber in the object, and on a linear combination of at least twomodel spectra, the model spectra representing wavelength-dependent basisfunctions of a radiation fluence or of a function of a radiation fluencein the object. Preferably, the radiation fluence is a normalizedradiation fluence and/or the model spectra represent the normalizedwavelength dependence of the radiation fluence.

The method for multispectral optoacoustic imaging of an object accordingto the invention comprises the following steps: irradiating the objectwith electromagnetic radiation at two or more different irradiationwavelengths, said electromagnetic radiation having a time-varyingintensity; detecting acoustic waves generated in the object uponirradiating the object with the electromagnetic radiation at thedifferent irradiation wavelengths; reconstructing images of the objectbased on the detected acoustic waves generated in the object at each ofthe irradiation wavelengths and determining a spatial distribution of atleast one first concentration value, which relates to a concentration ofat least one electromagnetic radiation absorber in the object, whereinthe determination of the spatial distribution of the at least one firstconcentration value is based on the reconstructed images at thedifferent irradiation wavelengths, on at least one wavelength-dependentextinction coefficient of the at least one electromagnetic radiationabsorber in the object, and on a linear combination of at least twomodel spectra, the model spectra representing wavelength-dependent basisfunctions of a radiation fluence or of a function of a radiation fluencein the object. Preferably, the radiation fluence is a normalizedradiation fluence and/or the model spectra represent the normalizedwavelength dependence of the radiation fluence.

The invention is based on the approach to determine a quantitativeimage, in particular a spatial distribution of a relative concentration,of at least one radiation absorber in the object considering i) two ormore optoacoustic images, which are reconstructed from acoustic wavesemanating from the object upon irradiating the object with transientelectromagnetic radiation at two or more different irradiationwavelengths, and ii) a linear combination of at least twowavelength-dependent model spectra that correspond to basis functions ofa, preferably normalized, wavelength-dependent radiation fluence withinthe object. Within the meaning of the invention a basis functionpreferably corresponds to an element of a particular basis for aspecific function space, whereby a linear combination of basis functionscan represent every continuous function in the specific function space(in analogy to the fact that every vector in a vector space can berepresented as a linear combination of basis vectors). Accordingly, inthe context of the invention, the term “model spectrum” and “basisfunction” is also referred to as fluence Eigenspectrum or Eigenspectrum.

According to the inventive approach, instead of using a model-basedoptical property estimation approach in which a model of lightpropagation in the object is used in order to extract radiation fluenceas well as absorption coefficients from images of an estimated absorbedenergy density in the object, a linear combination of a set ofwavelength-dependent model spectra is used, wherein the linearcombination is obtained by multiplying each model spectrum of the set ofmodel spectra by a wavelength-independent parameter and adding theresults. Preferably, the model spectra of a set of model spectra aredetermined in advance, e.g. by means of simulation of and/or measurementon a representative medium containing the at least one radiationabsorber, and are characteristic for the at least one radiationabsorber.

By means of the invention the following advantages, in particular overoptical property estimation approaches, are achieved: (i) A small-scaleoptimization problem with much smaller dimensionality can be defined.(ii) Application in a per-pixel basis is possible, meaning that full andaccurate image reconstruction of the whole illuminated volume is not anabsolute prerequisite. (iii) The application can be performed in aselected part of the whole image, whereby the convergence in anaccurately reconstructed area is not affected by possibly corruptedinformation in different parts of the image. (iv) The approach accordingto the invention relies on relative values of the absorbed energydensity (i.e. the wavelength dependence of the absorbed energy density:H(r,λ)/∥H(r)∥ whereby corresponds to the norm of the absorbed energydensity across all wavelengths at a position r) rather than absolutevalues of the absorbed energy density H(r,λ), avoiding furthercomplications like the need for perfect system calibration or estimatingthe Grüneisen coefficient of tissue, i.e. parameters which have beenshown to affect the performance of model-based optical propertyestimation approaches.

As a result, the invention allows for exact and reliable quantificationof a spatial distribution of at least one value relating to theconcentration of one or more chromophores in an object, in particularalso in tissue depths exceeding 3 mm.

In general, the invention relates to a device and a method forquantifying the distribution of one or more chromophores inside mediathat attenuate light using optoacoustic measurements. Although theinvention is particularly suited for quantifying tissue oxygenation, inparticular blood oxygen saturation sO₂, it can advantageously also beused in resolving and/or quantifying the relative concentration of anylight absorber distribution inside any medium containing at least onelight absorber.

According to a preferred embodiment, the processing unit is configuredto determine spatial distributions of at least two first concentrationvalues relating to concentrations of at least two electromagneticradiation absorbers, for example HbO₂ and HHb, in the object, whereinthe determination of the spatial distributions of the at least two firstconcentration values is based on the reconstructed images at thedifferent irradiation wavelengths, on the at least twowavelength-dependent extinction coefficients of the at least twochromophores, for example HbO₂ and HHb, present in the object, and alinear combination of the at least two model spectra, and to derive aspatial distribution of a second concentration value, in particular arelative concentration value like sO₂, from the determined spatialdistributions of the at least two first concentration values.

According to yet another preferred embodiment, the at least two modelspectra correspond to a minimized number of model spectra, which havebeen determined for the at least one electromagnetic radiation absorber,e.g. HbO₂ and HHb, assumed to be in the object. Preferably, theminimized number of model spectra is between 3 and 8, in particularbetween 3 and 5.

Moreover, it is preferred that the at least two wavelength-dependentmodel spectra have been determined by determining a set of spectralpatterns for different depths in a medium and for differentconcentrations of the at least one electromagnetic radiation absorber inthe medium, and by applying a statistical procedure on the set ofspectral patterns.

Preferably, the set of spectral patterns having been obtained fromanalytical solutions or from numerical solutions of the diffusionequation or the radiative transfer equation for different wavelengthsand/or from Monte Carlo simulations of light propagation for differentwavelengths and/or from a set of experimental measurements for differentwavelengths.

Preferably, the statistical procedure applied on the set of spectralpatterns comprises a principal component analysis (PCA) and/or a kernelprincipal component analysis (kernel PCA) and/or a linear discriminantanalysis (LDA).

According to another preferred embodiment, the linear combination of themodel spectra corresponds to an addition of a mean model spectrumΦ_(M)(λ) and a linear combination of a first number p of further modelspectra Φ_(i)(λ), wherein each of the further model spectra Φ_(i)(λ) ismultiplied by a model parameter m_(i): Φ_(M)(λ)+Σ^(p)_(i=1)m_(i)Φ_(i)(λ).

Preferably, the irradiation unit being configured to irradiate theobject with electromagnetic radiation at a second number w of differentirradiation wavelengths, wherein the second number w is larger than orequal to the sum of the first number p of further model spectra and anumber n_(c) of the first concentration values to be determined:w≧p+n_(c).

It is, moreover, preferred that the determination of the spatialdistribution of two first concentration values c′₁(r) and c′₂ (r), whichrelate to concentrations of two different electromagnetic radiationabsorbers in the object, comprises an inversion of a system ofnon-linear equations, wherein the inversion of the system of non-linearequations is performed simultaneously for a selected set of differentpositions r distributed in the image domain.

It is further preferred that the determination of the spatialdistribution of two first concentration values c′₁(r) and c′₂ (r), whichrelate to concentrations of two different electromagnetic radiationabsorbers in the object at a position r, comprises an inversion of asystem of non-linear equations given by

P(r,λ _(j))=(Φ_(M)(λ_(j))+Σ^(p) _(i=1) m _(i) ^(r)Φ_(i)(λ_(j)))(c′₁(r)ε₁(λ_(j))+c′ ₂(r)ε₂(λ_(j))), with j=1 . . . w,

P(r, λ_(j)) denoting the reconstructed image intensity of the object atan irradiation wavelength λ_(j) and position r, Φ_(M)(λ_(j))+Σ^(p)_(i=1)m_(i) ^(r)Φ_(i)(λ_(j)) denoting a linear combination of thewavelength-dependent model spectra Φ_(M)(λ_(j)) and Φ_(i)(λ_(j)), withi=1 . . . p, at an irradiation wavelength λ_(j) and position r, whereinm_(i) ^(r) denotes the values of the model parameters at a position r,and c′₁(r)ε₁(λ_(j))+c′₂(r)ε₂(λ_(j)) denotes a linear combination ofelectromagnetic radiation extinction coefficients ε₁(λ_(j)) andε₂(λ_(j)) of the two electromagnetic radiation absorbers at anirradiation wavelength λ_(j), wherein c′₁(r) and c′₂(r) denote the twofirst concentration values at a position r.

According to yet another preferred aspect of the invention, theinversion of the system of non-linear equations is performedsimultaneously for a number N_(grid) of different positions r_(n), withn=1 . . . N_(grid), corresponding to a grid in the image, and thevariation of one or more of the model parameters between neighborpositions is constrained.

Preferably, the system of non-linear equations is solved through theminimization of the difference between the reconstructed image intensityand the model (i.e. Σ_(j)∥P(r, λ_(j))−(Φ_(M)(λ_(j))+Σ^(p) _(i=1)m_(i)^(r)Φ_(i)(λ_(j)))(c′₁(r)ε₁(λ_(j))+c′₂(r)ε₂(λ_(j)))∥_(nm1) ^(pwr1), wherenrm1 denotes an appropriate norm and pwr1 an appropriate power) atmultiple positions corresponding to a grid in the image, and thesimultaneous minimization of the variation of model parameters betweenneighbor grid points. Preferably, this is achieved, e.g., through theminimization of a cost function f_(grid) that is a linear combination ofthe difference between the reconstructed image intensity and the modelat multiple positions corresponding to a grid in the image (i.e.Σ_(n)Σ_(j)∥P(r_(n), λ_(j))−(Φ_(M)(λj)+Σ^(p) _(i=1)m_(i)^(r(n))Φ_(i)(λ_(j)))(c′₁(r_(n))ε₁(λ_(j))+c′₂(r_(n))ε₂(λ_(j)))∥_(nrm1)^(pwr1)) and the simultaneous minimization of the variation of modelparameters corresponding to neighbor positions in the image (i.e. ∥m_(i)^(r(n1))−m_(i) ^(r(n2))∥_(nrm2) ^(pwr2) where r_(n1) and r_(n2)correspond to two neighbor positions in the grid and nrm2 denotes anappropriate norm and pwr2 an appropriate power). Preferably, theminimization function f_(grid) would take the following form:f_(grid)=Σ_(n) Σ_(j)∥P(r_(n), λ_(j))−(Φ_(M)(λ_(j))+Σ^(p) _(i=1)m_(i)^(r(n))Φ_(i)(λ_(j)))(c′₁(r_(n))ε₁(λ_(j))+c′₂(r_(n))ε₂(λ_(j)))∥_(nrm1)^(pwr1)+Σ_(n,m)λ_(n,m)∥m_(i) ^(r(n))−m_(i) ^(r(m))∥_(nrm2) ^(pwr2) whereλ_(n,m) are appropriate weight coefficients.

It is, moreover, preferred that in the inversion of the system ofnon-linear equations the values of the model parameters m_(i) ^(r) areconstrained to lie within a limited region (i.e. lim^(min) _(i)<m_(i)^(r)<lim^(max) _(i), where the values lim^(min) _(i) and lim^(max) _(i)are defined through simulations or experimental measurements), and/orthat the difference between the unknown values m_(i) ^(r) and somepredefined values {circumflex over (m)}_(i) ^(r) is also minimized alongwith the minimization of the cost function f_(grid). Preferably, this isperformed through the minimization of a second cost function g_(grid)that constitutes a linear combination of f_(grid) and a weighteddifference between the model parameters m_(i) ^(r) and the predefinedvalues {circumflex over (m)}_(i) ^(r) (also referred to as priors), i.e.g_(grid)=f_(grid)+w_(i) ^(r)∥m_(i) ^(r)−{circumflex over (m)}_(i)^(r)∥_(nrm3) ^(pwr3), where nrm3 denotes an appropriate norm, pwr3 anappropriate power and w_(i) ^(r) an appropriate weight coefficient.Preferably, the priors {circumflex over (m)}_(i) ^(r) are obtainedthrough simulations or experimental measurements.

The above and other elements, features, characteristics and advantagesof the present invention will become more apparent from the followingdetailed description of preferred embodiments with reference to theattached figures showing:

FIG. 1 an example of a device for multispectral optoacoustic imaging;

FIG. 2 an example of an Eigenspectrum model of wavelength-dependentoptical fluence: (a-d) Model spectra φ_(M)(λ), φ₁(λ), φ₂(λ) and φ₃(λ) atutilized wavelengths λ. (e) Mean Square Error (MSE) of the model on atraining dataset for different model dimensionalities. (f-h) Values ofparameters m₁, m₂ and m₃ on the training dataset as a function of depth(y axis) and tissue oxygenation (x axis);

FIG. 3 an example of an incorporation of constraints on the spatialcharacteristics of fluence parameters m₁, m₂ and m₃ in the case ofsimultaneous inversion on a grid of points. (a) A directed graph on thegrid of points enforces a constraint on the values of m₂ model parameterwith depth. (b) A non-directed graph on the grid of points penalizesgreat variations of the fluence parameters between neighbor points;

FIG. 4 an example of an sO₂ estimation error for all grid points usingthe Eigenspectrum method with simultaneous inversion on a grid of pointsand a non-negatively constrained linear unmixing approach.

FIG. 5 (a) an MSOT image (one wavelength presented) of a tumor-bearingnude mouse, wherein the points in the image correspond to the initialgrid selected for inversion. (b) Automatic refinement of the initialgrid for maximizing the SNR in the pixels used for inversion. (c) Map offluence parameter m₂ in the convex hull of the inversion grid. (d)Estimation of tissue blood sO₂ after fluence correction; and

FIG. 6 a schematic overview on an example of processes and algorithms ofa method for multispectral optoacoustic imaging.

FIG. 1 shows an example of a device for multispectral optoacousticimaging comprising an irradiation unit 10 which is configured toirradiate an object 1 with electromagnetic radiation at two or moredifferent irradiation wavelengths λ_(j) and a time-varying intensity.

For example, the irradiation unit 10 may comprise at least one lasersource which is configured to emit visible and/or near-infrared light,for example between 690 nm and 900 nm, preferably at 5 to 20 excitationwavelengths. Preferably, the emitted light is pulsed and/or amplitude-or frequency-modulated.

By irradiating the object 1 with transient, in particular pulsed and/ormodulated, electromagnetic radiation, acoustic waves are generated,which are detected by detection unit 20, which comprises one or moreultrasound detectors, e.g. one or more piezoelectric detectors.Preferably, the detection unit 20 comprises a two-dimensional orthree-dimensional array of ultrasound detector elements by which theacoustic waves generated in the object 1 can be detected and spatiallyresolved.

The object 1 can preferably be a biological tissue, like a part of ahuman or animal body. Depending on the particularities of the desiredapplication, pre-clinical or clinical applications, a holder (not shown)can be provided for holding and positioning the object 1, e.g. a tissuesample or a small animal, relative to the irradiation unit 10 and thedetection unit 20. Alternatively or additionally, the irradiation unit10 and detection unit 20 can be integrated in a handheld device whichcan be manually placed onto the object 1 and/or moved relative to theobject 1.

In the detector unit 20 the detected acoustic waves are converted intocorresponding detector signals which are forwarded to processing unit 30which is configured to reconstruct images of the object 1 based on thedetector signals that are obtained at each of the irradiationwavelengths λ_(j).

The processing unit 30 is further configured to determine a spatialdistribution of a relative concentration of at least one radiationabsorber contained in the object 1 by considering the reconstructedimages of the object 1 at different irradiation wavelengths λ_(j) andthe radiation fluence or a function of the radiation fluence in theobject 1, e.g. the normalized radiation fluence, wherein the radiationfluence or the function of the radiation fluence is based on a linearmodel of the radiation fluence. According to the linear radiationfluence model, at least two model spectra representingwavelength-dependent basis functions of the radiation fluence or afunction thereof are linearly combined.

Preferably, a display unit 40, e.g. a computer monitor, is provided fordisplaying the reconstructed images and/or the determined spatialdistribution of the relative concentration of the at least one radiationabsorber contained in the object 1.

In the following, preferred aspects of determining a spatialdistribution of the relative concentration of the at least one radiationabsorber contained in the object 1 will be elucidated in detail.

An algorithmic solution for obtaining quantitative measures of bloodoxygenation from multispectral optoacoustic images based on a newformulation of the quantification problem as a per-pixel non-linearspectral unmixing problem will be elucidated in more detail in section“PROBLEM FORMULATION” below.

Preferably, a simple linear model of wavelength-dependent basisfunctions is created for capturing the spectral perturbations of theoptical fluence in tissue using training data, which will be elucidatedin more detail in section “THE EIGENSPECTRUM MODEL OF OPTICAL FLUENCE”below.

Preferably, based on this model, a non-linear inversion scheme isformulated for obtaining quantitative blood sO₂ estimates from MSOTimages, which will be elucidated in section “INVERSION” below. Theproposed inversion scheme further exploits prior knowledge on structuraland spectral characteristics of the optical fluence for enhancinginversion accuracy.

Finally the application of the novel algorithmic solution on in vivomultispectral optoacoustic images will be demonstrated in section“APPLICATION ON TISSUE IMAGES” below.

Problem Formulation Optical Fluence and Acquiring Accurate sO₂ Estimates

With oxygenated (HbO₂) and deoxygenated (HHb) hemoglobin being the maintissue absorbers, optoacoustic images of tissues P(r,λ) can be relatedto the concentrations of HbO₂ and HHb through a non-linear relation,i.e.

$\begin{matrix}{{P\left( {r,\lambda} \right)} = {{C(r)}{{\Phi (r)}}_{2}{\frac{\Phi \left( {r,\lambda} \right)}{{{\Phi (r)}}_{2}} \cdot {\left( {{{c_{{HbO}\; 2}(r)} \cdot {ɛ_{{HbO}\; 2}(\lambda)}} + {{c_{HHb}(r)} \cdot {ɛ_{HHb}(\lambda)}}} \right).}}}} & (1)\end{matrix}$

In Eq. (1), C(r) is a space dependent scaling factor that includes theeffect of the Grüneisen coefficient (Γ(r)) and possibly additionalphysical phenomena (e.g. ultrasound attenuation) and system effects(e.g. the space and frequency dependent sensitivity of the ultrasoundsensors). ε_(HbO2)(λ) and ε_(HHb)(λ) are the wavelength dependent molarextinction coefficients of hemoglobin, while c_(HbO2)(r) and c_(HHb)(r)are the associated concentrations at a position r. ∥φ(r)∥₂ is the normof the optical fluence across all wavelengths at a position r, whileφ(r,λ)/∥φ(r)∥₂ is the normalized wavelength dependence of opticalfluence at a specific position, a term that will also be referred to asφ′(r,λ) in connection with the present invention.

The space-only dependent factors C(r) and ∥φ(r)∥₂ do not affect theestimation of sO₂ (Eq. (2)) which is calculated as a ratio once therelative concentrations of HbO₂ and HHb are known(c′_(HbO2)(r)=Const(r)·c_(HbO2)(r) and c′_(HHb)(r)=Const(r)·c_(HHb)(r),respectively).

$\begin{matrix}{{{sO}_{2}(r)} = {\frac{c_{{HbO}\; 2}^{\prime}(r)}{{c_{{HbO}\; 2}^{\prime}(r)} + {c_{HHb}^{\prime}(r)}}.}} & (2)\end{matrix}$

An estimate of the wavelength dependence of the optical fluence φ′(r,λ)is required for accurately extracting the relative values ofc′_(HbO2)(r) and c′_(HHb)(r). However, c′(r,λ) is in general an unknownfunction.

According to a preferred aspect of the invention, the wavelengthdependence of optical fluence φ′(r,λ) is not an arbitrary function ofwavelength, but it is assumed to only lie within a limited subspace ofspectral patterns within tissue. This aspect is based on the notion thatthe spectrum of optical fluence is a result of the cumulative lightabsorption by certain tissue absorbers, so the spectrum of opticalfluence will always be related to the spectra of these absorbers througha complex relation.

According to another preferred aspect, it is assumed that this subspaceof spectral patterns can be accurately described using a low-dimensionallinear model. In the section “THE EIGENSPECTRUM MODEL OF OPTICALFLUENCE” the creation of such a model that consists of a linearcombination of a small number of characteristic tissue optical fluencespectra, also referred to as Eigenspectra, will be described in detail.

Non-Linear Spectral Unmixing Formulation

According to yet another aspect, it is assumed that it is possible todefine a linear model that is based on a superposition of a meanspectrum φ_(M)(λ) and a number of optical fluence Eigenspectra φ_(i)(λ)to accurately characterize φ′(r,λ) at each arbitrary position r withintissue, then the non-linear mixture model for MSOT images at a positionr can be written as:

$\begin{matrix}{{P\left( {r,\lambda} \right)} = {\left( {{\Phi_{M}(\lambda)} + {\sum\limits_{i = 1}^{P}\; {m_{i}^{r}{\Phi_{i}(\lambda)}}}} \right).\left( {{{c_{{HbO}\; 2}^{\prime}(r)} \cdot {ɛ_{{HbO}\; 2}(\lambda)}} + {{c_{HHb}^{\prime}(r)} \cdot {ɛ_{HHb}(\lambda)}}} \right).}} & (3)\end{matrix}$

The solution for the P+2 unknowns (namely the P model parameters,m_(1 . . . P) ^(r) and the relative blood concentrations c′_(HbO2)(r)and c′_(HHb)(r)) can then be given using an appropriate non-linearinversion approach and at least P+2 excitation wavelengths. As statedbefore, the relative blood concentrations c′_(HbO2)(r) and c′_(HHb)(r)are proportional to the originals c_(HbO2)(r) and c_(HHb)(r) with regardto a scaling factor, however this fact does not affect the computationof sO₂, which is a ratiometric (relative) measurement.

In order to achieve the aforementioned non-linear spectral un-mixingformulation, a linear model for the optical fluence, in the form of(φ_(M)(λ), φ₁(λ), . . . , φ_(P)(λ)) is established. The model is able toaccurately describe the wavelength dependence of the optical fluenceφ′(r,λ), at any arbitrary position within tissue and independently ofthe surrounding optical properties and sO₂ values.

The Eigenspectrum Model of Optical Fluence

The wavelength-dependent optical fluence in tissue may take a number ofdifferent spectral patterns. The variation of optical fluence perwavelength is not random but follows a pattern typically determined bytissue depth, tissue scattering, blood sO₂ values and overall bloodconcentration. Preferably, by using Principal Component Analysis (PCA)on a set of different optical fluence spectral patterns, also referredto as training set, a linear model (φ_(M)(λ), φ₁(λ), . . . , φ_(P)(λ))consisting of a small number of fluence Eigenspectra φ_(i)(λ) isdetermined. PCA is used for offering a minimum square error property incapturing the spectral variability of the optical fluence spectralpatterns of the training set in a linear manner.

By testing the Eigenspectrum model in a set of optical fluence spectralpatterns produced in completely heterogeneous media with greatly varyingand randomly distributed optical properties and oxygenation values, itwas found that the Eigenspectrum model can accurately capture all suchoptical fluence spectral patterns with a very small mean square error(<2%), and it was confirmed that, independent of such conditions (i.e.optical property variation), the spectral patterns of optical fluencelie in a confined spectral subspace that can be accurately modeled bythe Eigenspectrum model.

In the following subsections, the light propagation models that arehereby used as well as an analytical description of the procedure forcreating the Eigenspectrum model from training data are elucidated indetail.

Light Propagation Models

Preferably, a light propagation model based on the diffusion equation isemployed, wherein the wavelength dependent optical fluence model(φ_(M)(λ), φ₁(λ), . . . , φ_(P)(λ)) is created using the solution of thediffusion equation for infinite homogeneous media, i.e.,

Φ_(z)=Φ₀exp(−μ_(eff) z),  (4)

where z denotes a distance from the source (tissue depth) and μ_(eff) isthe effective attenuation coefficient: μ_(eff)=−√{square root over(3μ_(α)(μ_(α)+μ_(s)′))}, μ_(s)′ being the reduced scattering coefficientand μ_(a) the absorption coefficient.

Alternatively, a finite element solution of the diffusion approximationfor arbitrary structure of optical properties may be used.

The Eigenspectrum Model

Preferably, for producing a training dataset to generate the opticalfluence Eigenspectrum model, simulations based on the analyticalsolution of the diffusion equation (Eq. (4)) are performed.

For example, a number of optical fluence spectral patterns φ_(z,ox)(λ)are computed for wavelengths λ=700 nm to 890 nm with a step size of 10nm, for depths from z=0 cm to z=1 cm with a step size of 0.143 mm, andfor uniform hemoglobin distribution at a constant oxygenation, varyingfrom ox=0% to ox=100% with a step of 5%. Constant optical properties ofμ_(α)=0.3 cm⁻¹ (at the isosbestic point of hemoglobin at 800 nm) andμ_(s)′=10 cm⁻¹ (constant for all wavelengths) are used. The so computed70×21 optical fluence functions φ_(z,ox)(λ) were normalized to disregardany scaling φ′_(z,ox)(λ)=φ_(z,ox)(λ)/∥φ_(z,ox)∥₂. The tissue depth of 1cm was selected for relating the application to small animaloptoacoustic imaging.

The computed vectors φ′_(z,ox) (λ) are used in the following as trainingdata in the context of PCA to derive a linear model of Eigenspectra (amean spectrum φ_(M)(λ) and three Eigenspectra φ₁(λ), φ₂(λ), φ₃(λ)) asexemplarily shown in FIG. 2(a)-(d).

PCA captures the dominant variations of all spectral patterns used inthe training dataset in the few first coefficients with a minimum meansquare error. A model dimension of 3 Eigenspectra was selected forproviding a rather simple model with a small associated error (see FIG.2(e)).

Regarding the parameters m₁, m₂ and m₃ of the model on the training dataset (see FIG. 2(f)-(h), respectively), it was surprisingly found thatthe values of m₂ can be greatly associated with tissue depth, while thevalues of m₁ can be greatly associated with the average levels ofbackground tissue oxygenation. This indicates that the secondEigenspectrum φ₂(λ) is mainly associated with the modifications ofoptical fluence due to depth and the average optical properties of thesurrounding tissue, while the first Eigenspectrum φ₁(λ) is associatedwith the “spectral shape” of optical fluence that relates to the averageoxygenation of the surrounding tissue. Moreover, the highest and lowestvalues of these parameters on the training dataset can offer anindication of their limits [lim^(min) _(i), lim^(max) _(i)] expected intissue.

A validation of the accuracy of the Eigenspectrum Model (φ_(M)(λ),φ₁(λ), . . . , φ_(P)(λ)) over simulation-based optical fluence spectracreated in arbitrary tissues yielded a very small model error,supporting the inventive approach according to which a simple linearmodel with only a few, e.g. three, Eigenspectra can capture the spectralvariability of φ′(r,λ) in complex tissue structures, independently ofthe optical properties, depth and blood oxygenation.

Similarly, the validity and accuracy of the Eigenspectrum model tocapture optical fluence spectral patterns present in tissue wasexperimentally investigated by means of measurements from small animalsin vivo and ex vivo, wherein the optical fluence in tissue was measuredby inserting a reference chromophore with well characterized spectrum indeep tissue. The result of this investigation yielded that theEigenspectrum model can capture realistically the spectral variabilityof optical fluence in tissue, even in the case of noisy experimentaldata obtained in vivo.

Inversion

Given the Eigenspectrum model for the optical fluence, the quantitativeestimation of blood sO₂ using MSOT reduces to solving Eq. (3) using anon-linear optimization approach. Inversion can be performedindependently for a single pixel in the image or simultaneously for agrid of pixels distributed appropriately in the image domain. In thesecond case additional constraints on the values of the Eigenspectrummodel parameters can be imposed, based on the spatial and spectralcharacteristics of optical fluence.

In the following, three different inversion approaches are presented,namely (1) a per-pixel inversion using local or global optimization, (2)a constrained optimization on a grid of points, where the spatialcharacteristics of the optical fluence are further modelled, and (3)another constrained optimization, wherein both spatial and spectralcharacteristics of the optical fluence are taken into account forenhancing inversion stability.

(1) Single Pixel Inversion

In the single pixel inversion case, the optimization problem for aspecific point r formulates as the minimization of the function in Eq.(5) over the values of m₁ ^(r), m₂ ^(r), m₃ ^(r), c′_(HbO2)(r) andc′_(HHb)(r), constrained according to Eq. (6).

$\begin{matrix}{f = {\begin{matrix}{{P\left( {r,\lambda} \right)} - {\left( {{\Phi_{M}(\lambda)} + {\sum\limits_{i = 1}^{P}\; {m_{i}^{r}{\Phi_{i}(\lambda)}}}} \right) \cdot}} \\\left( {{{c_{{HbO}\; 2}^{\prime}(r)} \cdot {ɛ_{{HbO}\; 2}(\lambda)}} + {{c_{HHb}^{\prime}(r)} \cdot {ɛ_{HHb}(\lambda)}}} \right)\end{matrix}}_{2}} & (5) \\{{\lim_{i}^{\min}{< m_{i}^{r} < \lim^{\max}}},{{c_{{HbO}\; 2}^{\prime}(x)} \geq 0},{{c_{HHb}^{\prime}(x)} \geq 0.}} & (6)\end{matrix}$

Using a constrained local optimization algorithm to solve the non-linearoptimization problem of Eq. (5) and (6) in the case of a simulated dataset it can be observed that the problem is non-convex and thus theinversion converges accurately for only a small percentage of pixels.Accordingly, the optimization problem is a so-called ill-posed problem.Preferably, in order to further improve the accuracy of the inversionand to achieve a unique solution, further constraints are imposed on thevalues of the fluence parameters m_(i) ^(r).

(2) Inversion on a Grid of Points Incorporating Spatial FluenceCharacteristics

Preferably, the spatial characteristics of optical fluence are furtherexploited for overcoming the ill-posed nature of the optimizationproblem. In contrary to tissue absorption which can vary arbitrarily,the optical fluence is bound to vary smoothly in space due to the natureof diffuse light propagation.

According to a preferred aspect of the Eigenspectrum model inversion,such a priori information is incorporated by attempting simultaneousinversion on a grid of points in the image domain, and penalizing largevariation of the fluence model parameters m₁ ^(r), m₂ ^(r) and m₃ ^(r)between neighbor pixels. This is preferably achieved throughincorporation of appropriate Lagrange multipliers λ_(i) to thenon-linear inversion function for constraining the variation of one orall model parameters. The values of the Lagrange multipliers can beoptimized using cross-validation on simulated data sets. Additionally,since the second Eigenspectrum is strongly associated with depth, the m₂model parameter can be constrained to take lower values as tissue depthincreases (in accordance with FIG. 2(g)).

FIG. 3 illustrates the incorporation of constraints on the spatialcharacteristics of fluence parameters m₁ ^(r), m₂ ^(r) and m₃ ^(r), inthe case of simultaneous inversion on a grid of points, wherein adirected graph on the grid of points enforces a constraint on the valuesof m₂ model parameter with depth (see FIG. 3(a)) and a non-directedgraph on the grid of points penalizes great variations of the fluenceparameters between neighbor points (see FIG. 3(b)).

Assuming a circular grid of P arcs and L radial lines (see FIG. 3) witha total of P×L points r_(p,l), and let the vector m_(i)=[m_(i) ^(r1,1),m_(i) ^(r1,2), . . . , m_(i) ^(r1,L), m_(i) ^(r2,1), . . . , m_(i)^(rp,l), . . . , m_(i) ^(rP,L)] correspond to the values of the opticalfluence parameter i (i=1 . . . 3) over all such points, the new inverseproblem is preferably defined as the minimization of function f_(grid)defined in Eq. (7) under the constraints defined in Eq. (8).

$\begin{matrix}{f_{grid} = {{\sum\limits_{i}\; {f\left( r_{i} \right)}} + {\lambda_{1}{{Wm}_{1}}_{2}} + {\lambda_{2}{{Wm}_{2}}_{2}} + {\lambda_{3}{{Wm}_{3}}_{2}}}} & (7) \\{{{\lim_{i}^{\min}{< m_{i}^{r_{k}} < \lim_{i}^{\max}}},{\forall i},{\forall k},{m_{2}^{r_{{p + 1},1}} < m_{2}^{r_{p,l}}},{m_{2}^{r_{{p + 1},{1 + 1}}} < m_{2}^{r_{p,l}}},{m_{2}^{r_{{p + 1},{1 - 1}}} < m_{2}^{r_{p,l}}}}{{{c_{{HbO}\; 2}\left( r_{k} \right)} \geq 0},{\forall k},{{c_{HHb}\left( r_{k} \right)} \geq 0},{\forall{k.}}}} & (8)\end{matrix}$

In Eq. (7), W is the weighted connectivity matrix corresponding to gridof points assumed. Each matrix element corresponds to a pair of gridpoints r_(p1,l1) r_(p2,l2) and is zero if the points are not directlyconnected or inverse proportional to their distance (w(r_(p1,l1)r_(p2,l2))∝1/∥r_(p1,l1)−r_(p2,l2)∥) if the points are connected.

Simultaneous inversion in a grid of points has surprisingly turned outto significantly enhance inversion stability.

(3) Inversion on a Grid of Points Incorporating Spatio-Spectral FluenceCharacteristics

For further enhancing the inversion stability, additional constraintscan be imposed to the values of fluence parameters based on priorexperimental measurements of the optical fluence in similar conditionsor based on alternative prior knowledge that can be derived directlyfrom the reconstructed multispectral optoacoustic images. Two possibleexamples of enforcing such constraints in the context of the gridinversion are described in the following:

Global Constraints on m₁

For instance, in certain cases of in vivo imaging, it can be preferablyassumed that the average tissue oxygenation is higher than 50%, and thusthe model parameter m₁ will obtain values lower or equal to zero asindicated by FIG. 2(f). This preferred assumption has been verified bymultiple controlled in vivo experiments (where the animal is breathing100% O₂ gas and imaging is performed in the brain or in the abdomen),where the measured fluence parameter m₁ was repeatedly found negative.Under this condition, and by constraining the values of m₁ ^(r) to benegative in all grid points excellent convergence rate can be achievedalso in the case of high levels superimposed noise. It is noted that noassumptions are made on the specific oxygenation values of the gridpoints, but only on the average tissue oxygenation.

Local Constraints on the Model Parameters

Alternatively, if no prior fluence measurements are available someconstraints for the model parameter m₁ can be extracted from theavailable data themselves. By performing a first linear unmixingapproach on the data a first estimation map of sO₂ levels can beobtained. It is noted that this sO₂ map is incrementally erroneous withtissue depth, however it can serve as a first approximation.

Using this sO₂ map and average tissue optical properties (e.g. μ_(α)=0.2cm⁻¹ at 800 nm and μ_(s)=10 cm⁻¹) the optical fluence within tissue canbe simulated and a first estimate on the values of one or all modelparameters m₁ ^(r), m₂ ^(r) and m₃ ^(r) (i.e. {circumflex over (m)}₁^(r), {circumflex over (m)}₂ ^(r), {circumflex over (m)}₃ ^(r)) can beextracted, where r corresponds to the positions of the grid points. Inthe following, the optimization problem of Eq. (7) to (8) is solved,with the only difference that the values of one or all of the modelparameters m₁ ^(r), m₂ ^(r) and m₃ ^(r) are constrained to lie in aregion close to the initial estimates {circumflex over (m)}₁ ^(r),{circumflex over (m)}₂ ^(r), {circumflex over (m)}₃ ^(r). It is notedthat these regions need to be incrementally larger with depth since indeep tissue the original estimates deviate usually from the true values.Such an optimization approach incorporating prior information for themodel parameters corresponding to each grid pixel can also be formulatedin the context of Bayesian inversion. By constraining the problem insuch a way the inversion stability is also greatly enhanced also incases of high random noise.

FIG. 4 relates to an example of an inversion of Eq. (7) and (8) on agrid of 64 points in the case of a simulated dataset and shows the errorin the estimation of local blood sO₂ for all grid points when using theinventive method (blue curve) and when using simple linearnon-negatively constrained unmixing (red curve). As apparent from thecomparison, the inventive method is associated with a very small errorin the estimation of blood oxygenation as compared to linear unmxingwhich is frequently used as an approximation.

Application on Tissue Images

In the following, the inventive Eigenspectrum method is exemplarilydemonstrated in the case of experimental MSOT data. A nude mouse bearingan orthotopic Mammarian 4T1 tumor was imaged in the tumor area in eightexcitation wavelengths from 690 nm to 900 nm.

FIG. 5(a) depicts the anatomical optoacoustic image acquired atwavelength 900 nm. The points on the image correspond to the initialgrid of pixels that were selected for the application of the non-linearinversion scheme described above. The grid was selected to cover thetumor area. The grid of pixels is in the following automaticallyadjusted to avoid using pixels of low image intensity (and thus low SNR)in the inversion.

After applying the non-linear inversion scheme on the refined grid ofpoints (see FIG. 5(b)), estimates for fluence parameters m₁ ^(r), m₂^(r) and m₃ ^(r) are obtained for each grid point r. In the following,the fluence parameters for each pixel in the convex hull of the grid isobtained using cubic interpolation, thus resulting in fluence maps forthe image sub-region.

FIG. 5(c) presents such a map for fluence parameter m₂. The wavelengthdependence of optical fluence is computed for each pixel within thesemaps as in φ′(r,λ)=φ_(M)(λ)+m₁ ^(r)φ₁(λ)+m₂ ^(r)φ₂(λ)+m₃ ^(r)φ₃(λ),where φ_(i)(λ) is the ith fluence Eigenspectrum.

Finally, a fluence-corrected MSOT image is obtained after diving theoriginal image P(r,λ) with the normalized wavelength dependent opticalfluence c′(r,λ) at each position r and wavelength λ:P^(corr)(r,λ)=P(r,λ)/φ′(r,λ). After fluence correction, non-negativelyconstrained linear spectral unmixing with the spectra of oxyanddeoxy-hemoglobin can result in an estimation of tissue bloodoxygenation.

FIG. 5(d) presents the estimated tissue oxygenation in the area of thetumor as computed using the Eigenspectrum method. Tissue oxygenation isoverlaid to the anatomical image with pseudocolor, represented in grayscale.

FIG. 6 shows an overview on an example of processes and algorithms of amethod for multispectral optoacoustic imaging. In present example,oxygenated (HbO₂) and deoxygenated (HHb) hemoglobin being the maintissue absorbers in the object 1.

In a first step, MSOT images are reconstructed at distinct irradiationwavelengths λ_(j) based on detected acoustic waves which are generatedin the object 1 (see FIG. 1) upon irradiating the object 1 withtransient electromagnetic radiation at the irradiation wavelengthsλ_(j).

In a second step, the reconstructed MSOT images undergo a pointselection process, wherein reliable pixels r_(i) in the image domain areselected in a semiautomatic way or by user annotation.

Optoacoustic image intensities P(r_(i),λ_(j)) at the selected pixelsr_(i) for all irradiation wavelengths λ_(j) are processed in a thirdstep, wherein fluence parameters m_(k) ^(ri) and concentration values,e.g. sO₂ values, are computed for each pixel r_(i) in the grid using theimage intensities P(r_(i),λ_(j)) and a pre-determined model of thefluence in the object, wherein the model is based on a linearcombination of model spectra φ_(k)(λ). This is performed by means of aninversion algorithm solving the following equation:

$\mspace{79mu} \begin{matrix}{\text{?}\left( {\text{?}{\sum\; \text{?}}} \right)\text{?}} \\{\text{?}\left( {\text{?}{\sum\; \text{?}}} \right)\text{?}} \\\ldots \\{\text{?}\left( {\text{?}{\sum\; \text{?}}} \right)\text{?}}\end{matrix}$ ?indicates text missing or illegible when filed

The calculated fluence parameters m_(k) ^(ri) (i=1 to Np, Np being thenumber of selected pixels, k=1 to K, K being the number of modelspectra) can be used to infer the wavelength dependence of the opticalfluence in the whole area corresponding to the convex hull of the gridby means of interpolation, i.e. the m_(k) values are interpolated fromm_(k) ^(ri) for all intermediate pixels in the convex hull of the gridof points.

Subsequently, the wavelength dependence of the optical fluence iscomputed for each pixel within the convex hull as in φ′(r,λ)=φ_(M)(λ)+m₁^(r)φ₁(λ)+m₂ ^(r) φ₂(λ)+m₃ ^(r)φ₃(λ)+ . . . .

A fluence-corrected MSOT image is obtained after dividing the originallyreconstructed MSOT image P(r,λ) with the normalized wavelength-dependentoptical fluence c′(r,λ) at each position r and wavelength λ:P^(corr)(r,λ)=P(r,λ)/φ′(r,λ).

After fluence correction, non-negatively constrained spectral unmixingwithin the spectra of oxy- and deoxy-hemoglobin result in thedetermination of a particularly accurate sO₂ map.

Preferably, the model of the fluence, in particular the model spectraφ_(k)(λ), is established independently, in particular in advance, of thecurrent measurements by means of a model creation algorithm, whereinbase functions of the wavelength dependence of the optical fluence inthe object are created for the excitation wavelengths λ_(j) employed.Preferably, in distinction to previous fluence modeling methods, onlythe wavelength dependence of optical fluence is modeled but not theactual fluence in absolute terms. Moreover, modeling is performed byapplying a statistical procedure on training data and not via a spatiallight propagation model that is fitted to the MSOT images.

SUMMARY

The invention relates to a novel solution for the problem ofquantitative multispectral optoacoustic imaging of tissue oxygenation.According to preferred aspects of the invention, a simple linear modelfor successfully capturing the spectral variability of optical fluencewith wavelength in tissue was created and validated. Using this model, anon-linear inversion scheme for correcting for the optical fluence andobtaining accurate quantitative estimates of tissue blood oxygenationwas formulated. Formulating the quantification problem of MSOT as aper-pixel non-linear unmixing problem instead of a model-based opticalproperty estimation approach allows for greater flexibility and easierapplicability in experimental images of tissue. The inventivemethodology does not assume full and accurate spatial information, anddoes not rely on absolute image values that are prone to error due amultitude of phenomena that are often not modelled in imagereconstruction (e.g. spatial sensitivity field of the ultrasound arrays,ultrasound attenuation and scattering, Grüneisen coefficient). In thisrespect the inventive method is readily applicable in in vivo MSOTimages. Furthermore, since the presented inverse problem is of lowdimensionality (typically 300-500 unknown parameters for the case ofgrid inversion), the computational cost is rather low. Specifically, therequired computational time for the inversion of a grid of 64 pointsusing an Intel Core i7-4820K can be less than 1 minute.

Based on a validation of the inventive approach, a rather simple3-dimensional model has turned out sufficient for accurately capturingthe spectral variability of optical fluence, independently of theoptical properties and the physiology of the surrounding tissue.

Basically, the absorption of water, lipids and melanin also present intissue may be ignored under an appropriate wavelength selection,assuming that the influence of hemoglobin is much more prevalent in thenear-infrared region. Moreover, the wavelength dependence of scatteringmay be ignored. However, according to preferred embodiments of theinvention, the influence of such parameters could also be incorporatedsimilarly for optimizing the Eigenspectrum model.

Despite the aforementioned simplifying assumptions, a model fit in thecase of in vivo and ex vivo experimental optoacoustic images of miceindicated a very promising match between theory and experimentalreality.

Because the inverse problem, in its simplest form (per pixel inversion),is non-convex and ill-posed, in a per pixel inversion case differentcombinations of optical fluence and local oxygenation may produce analmost identical result after multiplication. Thus, the nonlinearinversion may converge to a wrong solution unless additionalrestrictions are applied.

Therefore, it is particularly preferred to incorporate spatio-spectralfluence characteristics using available prior knowledge into a morecomplex inversion scheme. In this way, inversion stability was achievedeven in cases where the spectra are corrupted by random noise. Thepreferred inversion scheme provides accurate estimates of tissueoxygenation with far lower estimation error than the one of linearspectral un-mixing, that is typically used today in experimental studiesas an approximation.

What is claimed is:
 1. A device for multispectral optoacoustic imagingof an object, the device comprising a) an irradiation unit configured toirradiate the object with electromagnetic radiation at two or moredifferent irradiation wavelengths (λ), said electromagnetic radiationhaving a time-varying intensity, b) a detection unit configured todetect acoustic waves generated in the object upon irradiating theobject with the electromagnetic radiation at the different irradiationwavelengths (λ), and c) a processing unit configured to reconstructimages (P(r, λ)) of the object based on the detected acoustic wavesgenerated in the object at each of the irradiation wavelengths (λ) andto determine a spatial distribution of at least one first concentrationvalue (c′_(HbO2), c′_(HHb)), which relates to a concentration of atleast one electromagnetic radiation absorber (HbO₂, HHb) in the object,wherein the determination of the spatial distribution of the at leastone first concentration value (c′_(HbO2), c′_(HHb)) is based on thereconstructed images (P(r, λ)) at the different irradiation wavelengths(λ), at least one wavelength-dependent extinction coefficient(ε₁(λ_(j)), ε₂(λ_(j))) of the at least one electromagnetic radiationabsorber (HbO₂, HHb) in the object, and a linear combination of at leasttwo model spectra (Φ_(M)(λ), Φ_(i)(λ)), the model spectra (Φ_(M)(λ),Φ_(i)(λ)) representing wavelength-dependent basis functions of aradiation fluence (Φ(r,λ)) or of a function (Φ′(r,λ)) of a radiationfluence (Φ(r,λ)), in particular a normalized radiation fluence(Φ(r,λ)/∥Φ(r)∥), in the object.
 2. The device according to claim 1,wherein the processing unit is configured to determine spatialdistributions of at least two first concentration values (c′_(HbO2),c′_(HHb)) relating to concentrations of at least two electromagneticradiation absorbers (HbO₂, HHb) in the object, wherein the determinationof the spatial distributions of the at least two first concentrationvalues (c′_(HbO2), c′_(HHb)) is based on the reconstructed images (P(r,λ)) at the different irradiation wavelengths (λ), thewavelength-dependent extinction coefficients (ε₁(λ_(j)), ε₂(λ_(j))) ofthe at least two electromagnetic radiation absorbers (HbO₂, HHb) in theobject, and a linear combination of the at least two model spectra(Φ_(M)(λ), Φ_(i)(λ)), and to derive a spatial distribution of a secondconcentration value (sO₂(r)), in particular a relative concentrationvalue, from the determined spatial distributions of the at least twofirst concentration values (c′_(HbO2), c′_(HHb)).
 3. The deviceaccording to claim 1, wherein the at least two model spectra (Φ_(M)(λ),Φ_(i)(λ)) correspond to a minimized number of model spectra (Φ_(M)(λ),Φ_(i)(λ)), which have been determined for the at least oneelectromagnetic radiation absorber (HbO₂, HHb) assumed to be in theobject.
 4. The device according to claim 1, wherein the at least twowavelength-dependent model spectra (Φ_(M)(λ), Φ_(i)(λ)) have beendetermined by determining a set of spectral patterns (Φ_(z,ox)(λ)) fordifferent depths in a medium and for different concentrations of the atleast one electromagnetic radiation absorber (HbO₂, HHb) in the medium,and applying a statistical procedure on the set of spectral patterns(Φ_(z,ox)(λ)).
 5. The device according to claim 4, the set of spectralpatterns (Φ_(z,ox)(λ)) having been obtained from analytical solutions(Φ_(z,ox)(λ)=Φ₀ exp(−μ_(eff) z)) or numerical solutions of the diffusionequation or the radiative transfer equation for different wavelengths(λ) and/or Monte Carlo simulations of light propagation for differentwavelengths (λ) and/or a set of experimental measurements for differentwavelengths (λ).
 6. The device according to claim 4, the statisticalprocedure applied on the set of spectral patterns (Φ_(z,ox)(λ))comprising a principal component analysis (PCA) and/or a kernelprincipal component analysis (kernel PCA) and/or a linear discriminantanalysis (LDA).
 7. The device according to claim 1, the linearcombination of the model spectra (Φ_(M)(λ), Φ_(i)(λ)) corresponding toan addition of a mean model spectrum (Φ_(M)(λ)) and a linear combinationof a first number (p) of further model spectra (Φ_(i)(λ)), wherein eachof the further model spectra (Φ_(i)(λ)) is multiplied by a modelparameter (m_(i)): Φ_(M)(λ)+Σ^(p) _(i=1)m_(i) Φ_(i)(λ).
 8. The deviceaccording to claim 7, the irradiation unit being configured to irradiatethe object (1) with electromagnetic radiation at a second number (w) ofdifferent irradiation wavelengths (λ), wherein the second number (w) islarger than or equal to the sum of the first number (p) of further modelspectra (Φ_(i)(λ)) and the number (n_(c)) of the first concentrationvalues (c′_(HbO2), c′_(HHb)) to be determined: w≧p+n_(c).
 9. The deviceaccording to claim 1, wherein the determination of the spatialdistribution of two first concentration values (c′₁(r), c′₂(r)), whichrelate to concentrations of two different electromagnetic radiationabsorbers in the object, comprises an inversion of a system ofnon-linear equations, wherein the inversion of the system of non-linearequations is performed simultaneously for a selected set of differentpositions (r) distributed in the image domain.
 10. The device accordingto claim 1, wherein the determination of the spatial distribution of twofirst concentration values (c′₁, c′₂), which relate to concentrations oftwo different electromagnetic radiation absorbers in the object,comprises an inversion of a system of non-linear equations given byP(r,λ _(j))=(Φ_(M)(λ_(j))+Σ^(p) _(i=1) m _(i) ^(r)Φ_(i)(λ_(j)))(c′₁(r)ε₁(λ_(j))+c′ ₂(r)ε₂(λ_(j))), with j=1 . . . w, P(r, λ_(j)) denotingthe intensity of the reconstructed image of the object at a position rand an irradiation wavelength λ_(j), Φ_(M)(λ_(j))+Σ^(p) _(i=1)m_(i)^(r)Φ_(i)(λ_(j)) denoting a linear combination of thewavelength-dependent model spectra Φ_(M)(λ) and Φ_(i)(λ), with i=1 . . .p, at an irradiation wavelength λ_(j), wherein m_(i) ^(r) denotes modelparameters at a position r, and c′₁(r)ε₁(λ_(j))+c′₂(r)ε₂(λ_(j)) denotesa linear combination of electromagnetic radiation extinctioncoefficients ε₁(λ_(j)) and ε₂(λ_(j)) of the two electromagneticradiation absorbers at an irradiation wavelength λ_(j), wherein c′₁(r)and c′₂(r) denote the two first concentration values at a position r.11. A method for multispectral optoacoustic imaging of an objectcomprising the following steps: a) irradiating the object withelectromagnetic radiation at two or more different irradiationwavelengths (λ), said electromagnetic radiation having a time-varyingintensity, b) detecting acoustic waves generated in the object uponirradiating the object with the electromagnetic radiation at thedifferent irradiation wavelengths (λ), and c) reconstructing images(P(r, λ)) of the object based on the detected acoustic waves generatedin the object at each of the irradiation wavelengths (λ) and determininga spatial distribution of at least one first concentration value(c′_(HbO2), c′_(HHb)), which relates to a concentration of at least oneelectromagnetic radiation absorber (HbO₂, HHb) in the object, whereinthe determination of the spatial distribution of the at least one firstconcentration value (c′_(HbO2), c′_(HHb)) is based on the reconstructedimages (P(r, λ)) at the different irradiation wavelengths (λ), at leastone wavelength-dependent extinction coefficient (ε₁(λ_(j)), ε₂(λ_(j)))of the at least one electromagnetic radiation absorber (HbO₂, HHb) inthe object (1), and a linear combination of at least two model spectra(Φ_(M)(λ), Φ_(i)(λ)), the model spectra (Φ_(M)(λ), Φ_(i)(λ))representing wavelength-dependent basis functions of a radiation fluence(Φ(r,λ)) or of a function (Φ′(r,λ)) of a radiation fluence (Φ(r,λ)), inparticular a normalized radiation fluence (Φ(r,λ)/∥Φ(r)∥), in theobject.